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A model-independent analysis of the Fermi Large Area Telescope gamma-ray data 
from the Milky Way dwarf galaxies and halo to constrain dark matter scenarios 
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We implemented a novel technique to perform the collective spectral analysis of sets of multiple 
gamma-ray point sources using the data collected by the Large Area Telescope onboard the Fermi 
satellite. The energy spectra of the sources are reconstructed starting from the photon counts 
and without assuming any spectral model for both the sources and the background. In case of 
faint sources, upper limits on their fluxes are evaluated with a Bayesian approach. This analysis 
technique is very useful when several sources with similar spectral features are studied, such as 
sources of gamma rays from annihilation of dark matter particles. We present the results obtained 
by applying this analysis to a sample of dwarf spheroidal galaxies and to the Milky Way dark matter 
halo. The analysis of dwarf spheroidal galaxies yields upper limits on the product of the dark matter 
pair annihilation cross section and the relative velocity of annihilating particles that are well below 
those predicted by the canonical thermal relic scenario in a mass range from a few GeV to a few 
tens of GeV for some annihilation channels. 
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I. INTRODUCTION 

Milky Way dwarf spheroidal (dSph) galaxies are 
candidate targets for dark matter (DM) studies through 
annihilation signatures. This is because their mass-to- 
light ratio is predicted to be of the order of 10 — 10 3 [l|,0, 
implying that they could be largely DM dominated. 
Moreover, since no significant gamma-ray emission of 
astrophysical origin is expected (these systems host few 
stars and no hot gas) , the detection of a gamma-ray signal 
could provide a clean DM signature. 

The Milky Way halo is another promising candidate 
for DM searches. An approach to search for DM 
emission from annihilation in the Galactic halo is to 
study the gamma-ray flux from sky positions distant 
from known astrophysical gamma-ray sources. The 
diffuse emission from unresolved sources and from the 
interaction of charged particles with the interstellar 
medium constitutes a background for this approach. 

Weakly Interacting Massive Particles (WIMPs) have 
long been considered as well-motivated candidates for 
DM that could contribute to the 80% of the non-baryonic 
mass density in the universe [3j . 

At a given energy E, the differential gamma-ray flux 
$ 7 (£', Ail) (in units of photons cm -2 s _1 GeV -1 ) from 
WIMP annihilation in a region covering a solid angle Ail 
and centered on a DM source, can be factorized as Pf: 



$ 7 (£, Afl) = J(Afi) x $ pp (£) (1) 

where J (Aft) (in units of GeV 2 cm -5 sr) is the 
"astrophysical factor" (hereafter, J-factor), i.e., the line 
of sight (l.o.s.) integral of the DM density squared in the 
direction of observation over the solid angle Afl: 
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The term <& PP (E) (in units of GeV -3 cm 3 s -1 sr -1 ) 
is the "particle physics factor", that encodes the particle 
physics properties of the DM, and for a given WIMP 
mass m x is given by: 
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where (av) is the WIMP pair annihilation cross section 
times the relative velocity of the two annihilating 
particles, while Bf and Nf(E,m x ) are respectively the 
branching ratio and the differential photon spectrum of 
each pair annihilation final state /. 

We note that the particle physics factor in Eq. [1] is 
independent of the spatial distribution of the DM, and 
hence independent of the particular DM source under 
investigation. Eq. [T]can be rewritten as: 
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showing that the ratio between the photon flux and the 
J-factor is expected to be independent on the source if 
the DM annihilation mechanism is the same for all the 
sources. Starting from a measurement of the gamma- 
ray flux from a candidate DM source, if the J-factor is 
known, Eq. |4] allows us to obtain a measurement of the 
particle physics factor. If the kinematic terms of the 
summation in Eq. [3] are known, this measurement will 
yield an estimate of (av) as a function of m x . Moreover, 
since <& PP (E) is independent of the source, the results 
from individual sources can be combined, thus improving 
the measurement. 

Recently two analysis approaches were developed to 
constrain DM models using the Fermi LAT data [a, Q . 
In Ref. [5|, a binned Poisson likelihood fit was used to 
fit both the spatial and the spectral information for the 
reconstructed photon events collected in a sky region 
around the target source. The data from 10 dSphs were 
also combined using a joint likelihood analysis that takes 
into account the uncertainties on the J-factors. The 
upper limits on (av) were evaluated implementing an 
approach based on a profile likelihood function, that 
incorporates the nuisance parameters. In Ref. [5j a 
two-year photon data sample was analyzed with the 
P6_V3_DIFFUSE Instrument Response Functions (IRFs) 
in the energy range from 200 MeV to 100 GeV. In Ref. [6[ 
a three-year photon data sample was analyzed with the 
P7SOURCE_V6 IRFs in the energy range from 1 to 
100 GeV. Photons from a sky region with an angular 
radius of 0.5° from each dSph were selected and the 
background was evaluated by sampling positions within 
an angular distance of 10° from each dSph and counting 
the number of events in a cone of 0.5° angular radius. 
The upper limits were evaluated using a fully frequentist 
approach that takes in account the different J-factors of 
each source. The authors also took the uncertainties on 
the J-factor into account with a semi-Bayesian approach 
as well. The results of these two analyses were used to set 
upper limits on the annihilation cross section below the 
canonical value for the thermal relic WIMP scenario of 
3 x 10_- 26 cm 3 s" 1 [1 0] up to masses of about 30 GeV for 
the bb and t + t~ channels. This limit could represent a 
serious challenge to the conventional WIMP dark matter 
hypothesis. 

In this work we present the results obtained with 
a model-independent data analysis method 8] applied 
to DM searches. This method can be applied to the 
analysis of individual sources (Sect. IIII Aft as well as 
to the combined analysis of multiple sources (Sect.s 
IIIIBI and IIIICp . and does not introduce degrees of 
freedom in the calculation of confidence intervals on 
the parameters in the DM model. The first step of 
the analysis is the evaluation of upper limits on the 
possible gamma-ray signal events. This calculation is 
performed by properly choosing, for each source, a signal 
and a background region (see Sect.s [ill] and IIV|) and 
following a Bayesian approach to evaluate upper limits 
on the signal counts. In this way no models are required 



for the source and for the background. Moreover, the 
effects of systematic uncertainties can be easily taken 
into account by integrating over a nuisance parameter 
(either the J-factor or the effective area) the posterior 
probability distributions (Sect. IIII D[) . Finally, the upper 
limits on the photon counts can then be converted into 
upper limits on & PP (E), and consequently on (av), once 
a DM model has been implemented. In the present 
analysis we used a sample of gamma-ray data collected 
by the Fermi LAT during its first 3 years of operation 
in survey mode. The data were analyzed using the most 
recent LAT IRFs (P7SOURCE_V6 and P7CLEAN.V6). 
Candidate photons converting in both the front and back 
part of the instrument in the energy range from 562 MeV 
to 562 GeV were used for the analysis. Upper limits on 
(av) as a function of m x were obtained from the analysis 
of individual dSph galaxies (Sect. IIV|) and from their 
combined analysis, as well as from the analysis of the 
Milky Way Halo (Sect. EJ). 



II. THE INSTRUMENT AND THE DATA 

The LAT is a pair-conversion gamma-ray telescope 
designed to measure gamma rays in the energy range 
from 20 MeV to more than 300 GeV. In this paper a 
brief description of the LAT is given, while full details 
can be found in (l2| |. 

The LAT is composed of a 4 x 4 array of 16 
identical towers designed to convert incident gamma- 
rays into e + e~ pairs, and to determine their arrival 
directions and energies. Each tower hosts a tracker 
module and a calorimeter module. Each tracker module 
consists of 18 x-y planes of silicon-strip detectors, 
interleaved with tungsten converter foils, for a total 
on-axis thickness equivalent to 1.5 radiation lengths 
(r.L). Each calorimeter module, 8.6 r.l. on-axis thick, 
hosts 96 CsI(Tl) crystals, hodoscopically arranged in 
8 perpendicular layers. The instrument is surrounded 
by a segmented anti-coincidence detector that tags the 
majority of the charged-particle background. 

A sample of gamma-ray data collected by the Fermi 
LAT during its first three years of operation in 
survey mode was used for this analysis, overlapping 
substantially with the data used for the second LAT 
source catalog [13J . The Pass7 IRFs [14| event selection 
cuts (for SOURCE and CLEAN event classes), with 
candidate photons converting in both the front and 
back parts of the instrument, were used. To avoid 
contamination from the bright limb of the Earth, 
data taken during any time period when the angular 
separation of a cone of 10° angular radius centered on 
the source direction with respect to the Zenith direction 
exceeded 105° were discarded, as well as data taken 
during any time period when the LAT rocked to an angle 
exceeding 52°. The data taken during time periods when 
the source was observed with an off-axis angle larger than 
66.4° were also discarded. 



We performed the spectral analysis using the inter- 
nal LAT Collaboration software package FermiUnfold- 
ing [SHllI , which enables gamma-ray spectra to be recon- 
structed without assuming any model for the sources or 
the background. The data analysis was performed select- 
ing gamma rays with energies from 562 MeV to 562 GeV. 
The energy interval was divided into 12 bins, equally 
spaced on a logarithmic scale (4 bins per decade). We 
emphasize that, to take energy dispersion into account, 
in the unfolding approach there is a distinction between 
the observed photon energies and the true ones. The re- 
lationship between observed and true energy is expressed 
in terms of a smearing matrix, which represents the IRF 
and is evaluated by means of a full Monte Carlo simula- 
tion m. 



livetimes of the signal and of the background regions. 
However, since the data selection cuts illustrated in fJTT] 
are performed on a cone of 10° angular radius centered on 
the source, and since the outer radius of the background 
annulus used for the present analyses is always less than 
10° (see PV[) . the livetime ratio T s /Tb is always equal to 
1. 

The posterior probability density function (PDF) of 
the signal counts s was calculated assuming a uniform 
prior PDF for both s and b and is given by [15j: 



p(s\n,m) = \ OfcS e 



(6) 



fc=0 



with the coefficients a^ defined as: 



III. ANALYSIS METHODS 
A. Study of individual sources 

For each individual source a signal region and a 
background region were defined. The signal region, in 
which gamma rays emitted from the source are expected, 
was defined as a cone of a given angular radius, centered 
on the nominal position of the source. On the other hand, 
the background region was usually defined as an annulus 
centered on the source position and external to the 
signal region. To rule out possible contaminants in the 
background evaluation, when defining the background 
regions all the sources in the 2FGL catalog [13J were 
masked. The values of the angular radii adopted in this 
analysis to define the signal and background regions, as 
well as for masking the 2FGL catalog sources, are given 
in §|TV]andin § [V] 

Since the possible gamma-ray signal is expected to be 
faint, in each energy energy bin we set upper limits on 
the signal counts. The evaluation of the upper limits was 
performed following the Bayesian approach illustrated 
in Ref. (la|. Following the notation of Ref. [15[, we 
indicate with n and m the number of photons detected 
in a given energy bin in the signal and background 
regions, respectively (in the following, to keep the 
notation simple, we will suppress the energy dependence 
of these variables). We assume that the probabilities 
of measuring the pair (n,m) are both Poissonian with 
expectation values s + cb and b, respectively, where s is 
the expectation value of the signal counts (in the signal 
region), b is the expectation value of the background 
counts (in the background region) and c is defined as: 



aK 



(5) 



where Af2 Sj {, are the solid angles of the signal and of 
the background regions respectively. In principle the 
definition of c in the previous equation should include 
the livetime ratio T s /T{,, where T St b are respectively the 
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where J\f is a normalization constant. 

In case of the absence of a background (c = 0, m = 0), 
the posterior PDF on the signal reduces to (l5| : 



p(s\n) 



r(n + l)' 



(8) 



The upper limit on the signal counts s u at the 
confidence level (or credibility level, CL) 1 — a was 
evaluated by numerically solving the integral equation: 



p(s\n, m)ds = 1 



(9) 



The upper limits on the signal counts were finally 
converted into upper limits on the flux by means of the 
unfolding procedure described in [9l4ll|. The smearing 
matrix associated with each sky direction was built by 
taking into account the pointing history recorded by 
the LAT Q and was evaluated from the Monte Carlo 
simulation of the LAT. 

The measured upper limits on the flux were then 
converted into upper limits on (<rv). From Eqs. [3] and 0] 
it follows that: 
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(10) 



For each value of m x the conversion of the limits on the 
gamma-ray flux into limits on (av) was performed by 
requiring that the flux predicted from the model must 
not exceed the measured photon flux in any energy bin. 
The expected gamma-ray flux from the DM annihilation 
channels was evaluated as a function of energy using 
the DMFIT package [lj| based on DarkSUSY [nj], as 
implemented in the LAT Science Tools [18|. For large 



DM masses (around or above ITeV), the radiation of 
soft elcctroweak bosons leads to additional gamma rays in 
the ener gy r ange of relevance for the present analysis (see 
e.g. [ill [20| ) . This emission mechanism is not included 
in the DMFIT package. Therefore the present analysis 
provides conservative upper limits on (crv}. 



B. Stacking analysis 

According to Eq. [H the particle physics factor is 
independent of the source under investigation. This 
feature suggests the possibility of combining the data 
from all individual sources in order to improve the 
constraints on the DM models. 

Once the individual sky directions were analyzed, a 
stacking analysis was performed. In this case the counts 
from the signal and background regions corresponding 
to each source were added, and the upper limits on the 
signal were evaluated following the same procedure as for 
individual sources. 

In order to implement the same analysis procedure as 
for individual sources, in the stacking analysis the ratio 
between the signal and background regions was defined 
as: 



to stacking the events from each sky direction on top of 
one another and then analyzing the resulting image (as 
an example see the last panels in Figs. [5] and [3] for the 
case of the dSph analysis) . 

Since the sources may be seen by the instrument 
with different exposures, the J- factor value used in 
the stacking analysis was defined as the average of 
the J- factors of individual sources, each one weighted 
with its exposure in the whole energy interval under 
investigation. In principle, different J-factors should be 
determined for each energy bin, with each one evaluated 
taking into account the exposures in the corresponding 
bin. We performed this calculation in the case of the 
dSph galaxies, and the differences between the J-factors 
evaluated using the exposures in individual energy bins 
with respect to the J-factor evaluated using the overall 
exposure were less than 1%. Since these differences are 
small, we decided to use the same J-factor for the whole 
energy interval, evaluated using the overall exposure. 
In this way we also avoided introducing an energy 
dependence of the J-factor that may seem unphysical 
since, according to Eq. [5J the J-factor is determined only 
by the DM density profile. 

The upper limits on (<rv) were evaluated in the same 
way as for individual sources. 
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where Afigj^j and T a i t tn are respectively the solid angles 
and the livetimes of the signal and background regions of 
the i-th source (T s i = Tu according to the discussion in 

MM- 

We note that a more detailed statistical analysis (see 
the discussion in Appendix |Aj shows that the coefficient 
c should be defined as: 
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where Cj is the coefficient defined in Eq. [S] for the i-th 
source and mi are the counts in the background region 
of the z-th source. 

If the coefficient c is defined as in Eq. [TTJ its value 
will depend only on the extensions (solid angles) of the 
signal and background regions and on the livetimes of 
the stacked sources. On the other hand, if c is defined as 
in Eq. [TJ1 its value will also depend on the data (counts 
in the individual background regions). In the stacking 
analysis of the dSph galaxies we evaluated the coefficients 
c using both the definitions in Eq. QT] and Eq. [TJJ We 
found that the differences in the values of c obtained 
implementing the two different definitions were negligible 
in all the energy bins. 

The different exposures of the individual sources were 
also taken into account in the evaluation of the smearing 
matrix [9| , which was performed by stacking the pointing 
histories of all the sources. This procedure is equivalent 



C. Composite analysis 

In the previous analysis the events from all the sources 
were stacked. This is equivalent to considering the set of 
all the sources as a single source with a J-factor given 
by the average value weighted with the exposures of all 
the sources. In the stacking method all the sources are 
treated in the same way, and photons from a source 
with a small J-factor are considered as likely to originate 
from DM as photons from a source with a large J-factor. 
However, in the absence of a clear gamma-ray signal, 
i.e., if the counts in the signal region n are compatible 
with the expected background (n ~ c m), a source with a 
higher J-factor will yield a lower upper limit with respect 
to a source with a lower J-factor. The "DM sensitivity" 
of each source is therefore determined by its J-factor. 

To account for the different sensitivities of each 
observation we developed a composite analysis approach 
that combines the results from all the sources taking into 
account the individual J-factors. Unlike the approaches 
discussed in sections flll Al and lHI B[ for simplicity, in this 
approach we do not treat the energy resolution. Since the 
68% containment of the energy resolution of the LAT 
in the energy range chosen for the present analysis is 
less than 15% [lj], we expect that neglecting the energy 
dispersion in the evaluation of the flux could yield a 
similar uncertainty. Indeed we verified, in the case of 
dSph galaxies, that the differences between the fluxes 
reconstructed either neglecting or taking into account the 
energy dispersion are of the order of a few percent in the 
whole energy range of the analysis. 



Indicating with Si the expected signal counts from the 
i-th source in the energy interval [E, E + AE] (again, 
for simplicity, we will suppress the energy dependence 
of these variables), it is possible to define the random 
variable u as: 



u = riiSi 
with the factor r\i defined as: 



(13) 



Expanding the calculations in the previous equation, the 
final expression of the likelihood function is given by: 
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where r\ and the set of coefficients //. are defined as 
follows: 



1 



m = 
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where J{ is the J- factor of the z-th source and £i (E) is its 
exposure in the energy bin [E, E + AE], which is given 
by: 



(15) 



£i(E)= / dtf LT (t)ME,t) 



where Ai(E,t) is the effective area and fLr(t) is the 
livetime fraction. The dependence on t in Eq. ll5l indicates 
that the aspect angles (off-axis and azimuthal angles in 
the instrument frame) corresponding to the given sky 
direction (source) are changing with time. 

Since u is equal to (1/J)<fr 7 (i?, Afi), and hence to 
$ {E), it is independent of the particular source under 
investigation. 

A set of PDFs for the random variable u can 
be evaluated starting from the data of each source. 
Indicating with n, and m, the counts in the signal and 
background regions of the i-th source, the PDF for s, 
is given by Eq. [51 which can be rewritten explicitly 
indicating the source index as: 



P 
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with the coefficients a^. defined as in Eq. [7] 

The i-th PDF for the variable u can be derived from 
Eq. 1161 and is given by: 



p l {u\n i ,m l ) = e u/vi 2J b ik i 
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with the coefficients 6,i- defined as 
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To combine a set of N sources we build the likelihood 
function: 



N 

C(u\n 1 ,mi;n 2 , m 2 ; . . . ; n N , m N ) = J^Pi(u|ni, m»). 
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In the summation of Eq. [20] the maximum value of k 
yielding a non-zero coefficient fk is n„ laa; = ni + n-i + 
.. . + n N . 

In case of the absence of a background the expression 
of the likelihood function becomes simpler. Starting from 
Eq. [51 it is straightforward to show that the expression 
of the likelihood function is given by: 
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where 



JV 



(24) 



The likelihood function obtained from Eq. [2D] (or 
Eq. I23p is not normalized because, having assumed 
that u — Si/rji is independent of the source 
under investigation, the measurements (rij,TOj) are not 
independent of each other. To ensure normalization, the 
function C(u) must be multiplied by a constant A, which 
in the general case of Eq. [2UJ is given by: 



A = 
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Once the likelihood function is normalized, the upper 
limit u* at the CL 1 — a can be evaluated by numerically 
solving the equation: 



AC{u)du = 1 



a. 



(26) 



D. Systematic uncertainties 

Systematic uncertainties on the J-factor as well as 
on the effective area can be taken into account in 
the above procedures introducing a nuisance parameter 
in the definition of the random variable u. In the 
following we will illustrate the calculations to take into 
account the systematic uncertainties on the J-factors; the 
mathematical formalism used in the calculations to take 
into account systematic uncertainties on the effective 
area is similar so we do not present it here. 

Similar to the approach in §111 CI it is possible to define 
the random variable u as: 



Ji 



where the the factor pi is defined as: 
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(27) 



(28) 



Unlike in mil CI in this case the dependence of u on J; is 
written explicitly in order to take fluctuations in J^ into 
account. 

The posterior PDF for u can be obtained starting from 
the joint PDF p,(sj, Ji) for Si and Ji as: 



Pi(u) = — Ji Pi (Jiu/pi, Ji) dJi. (29) 

Pi J 

Since Si and Ji are independent random variables, their 
joint PDF can be factorized, and the previous equation 
rewritten as: 



area of the instrument. In this case the effective area 
is treated as a uniformly distributed random variable, 
while the J-factor is assumed to be known. 

For the stacking analysis, the same procedure was 
implemented as for the analysis of individual sources. In 
this case the PDF for u was obtained starting from the 
cumulative signal and background counts, as in §111 Bl 

In the case of the composite analysis, the likelihood 
function was built by multiplying all the individual PDFs 
Pi(u) in Eq. 1311 and then the upper limits on u were 
evaluated as discussed in 3111 CI 



IV. ANALYSIS OF THE DWARF SPHEROIDAL 
GALAXIES 

This analysis was performed using P7SOURCE_V6 
class events. For each source, events within a cone of 
10° angular radius centered on the nominal position of 
the source were selected. Again we note that, because 
of the selection cuts described in §TH all sky directions 
within 10° from the source will have the same live time. 

The positions of the dSph galaxies considered in the 
present analysis and the corresponding values of their 
J-factors are reported in Tab. U which is taken from 
Ref. [5J]. These dSph galaxies are not included in the 
second catalog of the Fermi LAT [l3|, i.e., they are not 
detected in the gamma-ray energy band above 100 MeV. 
In our analysis we assumed that the J-factor distribution 
for each dSph is well described by a log-normal function 
(see Ref. [a] for more details), with average value and 
standard deviation of log 10 J reported in Tab.Q] The half- 
light radii of the dSph galaxies used to compute the J- 
factors are less than or close to 0.5°. The average values 
(Ji) were calculated from the log-normal distributions as: 



Pi(u) = — Ji pi{Jiu/pi) Pi(Ji) dJi. (30) 

Pi J 

where the PDF pi(si) is given by Eq. [16] 

To make the calculation simpler, for the J-factors 
a uniform PDF in the range [Jn, Ji2\ is assumed, i.e. 
Pi(Ji) — 1/AJj. Introducing these PDFs in the previous 
equation, the posterior PDF for u is given by: 
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The upper limit u* at the CL 1 — a is evaluated by 
numerically solving the integral equation: 



Pi{u)du = 1 



(32) 



where Pi(u) is given by Eq. 1311 

A similar approach is implemented to evaluate the 
effects of the systematic uncertainties on the effective 
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where pi and of are the average value and the variance 
of the distributions of In Ji , which can be calculated 
multiplying the values reported in Tab. Uby In 10. 

When evaluating the effects of systematic uncertainties 
on the J-factor, we followed the procedure discussed 
in mil Dl The fluctuations of Ji were described by a 
uniform PDF in the interval [Jn, J 12], with Jn and J^ 
corresponding to the 16% and 84% quantiles of the log- 
normal distribution. 

In the case of the stacking analysis, as discussed 
in 3111 Bl the J-factor was evaluated as the weighted 
average of the J-factors of individual dSph galaxies 
with the exposures. The distribution of the J-factors 
of the stacked sources was built sampling a large 
set of events (10 6 ) from the J-factor distributions of 
individual sources, and is shown in Fig. [1] The 
average value of the J-factor for the stacked sources 
is (J) = 1.75 • 10 19 GeV 2 cm~ 5 sr. As in the case of 
individual sources, to study the systematic uncertainties 



Name 


Galactic 


Galactic 


log 10 (J) 




longitude 


latitude 


( GeV 2 cm- 5 sr) 


Bootes I 


358.08° 


69.62° 


17.7 ±0.34 


Carina 


260.11° 


-22.2° 


18.0 ±0.13 


Coma Berenices 


241.9° 


83.6° 


19.0 ±0.37 


Draco 


86.37° 


34.72° 


18.8 ±0.13 


Fornax 


237.10° 


-65.7° 


17.7 ±0.23 


Sculptor 


287.15° 


-83.16° 


18.4 ±0.13 


Segue I 


220.48° 


50.42° 


19.6 ±0.53 


Sextans 


243.4° 


42.2° 


17.8 ±0.23 


Ursa Major II 


152.46° 


37.44° 


19.6 ± 0.40 


Ursa Minor 


104.95° 


44.80° 


18.5 ±0.18 



TABLE I: List of the dSph galaxies used in this analysis. 
The J-factors are assumed to be distributed according to a 
log-normal distribution with (log 10 J) and o"io glo j given here. 
The half-light radii of the dSph galaxies used to compute the 
J-factors are less than or close to 0.5° [5J. 



on J we used a uniform PDF in the interval [ J\ , J^\ , with 
J\ and Ji corresponding to the 16% and 84% quantilcs 
of the resulting J-factor distribution. 

For each dSph galaxy the signal region was defined 
as a cone of angular radius A6 = 0.5° centered on the 
source position. The value of A8 is the same as the one 
used to evaluate the J-factor, and is consistent with the 
Point Spread Function (PSF) of the instrument, the 68% 
containment radius of which is smaller than 1° in the 
energy range above 1 GeV [lj] . The background region 
was defined as an annulus centered on the source position, 
with an inner radius of 5° and an outer radius of 6°. In 
order to prevent contamination of the background sample 
from photons emitted by other astrophysical sources, all 
the sources in the 2FGL Catalog [13| were masked. Using 
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FIG. 1: Distribution of 10 6 J-factor values for the stacked 
sources. Each realization is obtained by sampling the 10 log- 
normal distributions of the J-factors of individual sources and 
evaluating the average value weighted by the exposures. The 
red lines correspond to the 16% and 84% quantiles of the 
distribution. 



the HEALPix |2l| pixelization scheme with A^ide = 256, 
the sky was divided into 786432 equal area pixels, each 
covering a solid angle of 1.6 • 10~ 5 sr. The background 
region was composed of all the pixels in the annulus, 
excluding those at an angular distance less than 3° from 
any point source and those at an angular distance less 
than 3° plus the angular size of the semi-major axis from 
any extended source. The solid angle of the background 
region, Ailti, was then evaluated by adding the solid 
angles corresponding to all the unmasked pixels in the 
annulus. 

Fig. [2] shows the photon count maps with energy 
greater than 562 MeV for the 10 dSph galaxies considered 
in the present analysis. A qualitative inspection of the 
count maps shows no evidence of a gamma-ray signal 
from any dSph galaxy. On the other hand, from Fig. [2j 
bright gamma-ray sources close to some dSph galaxies are 
evident. However, as mentioned above, these sources are 
not considered when evaluating the background because 
of the masking procedure. 

Photons emitted by possible gamma-ray point sources 
lying close to a dSph galaxy might be detected in the 
signal region. These photons will not be accounted for in 
the background and, therefore, they might be confused 
with a DM annihilation signal. The result is that the 
upper limits on the DM signal will be higher and so, in 
this sense, this analysis is conservative. 

Fig. [3] shows the signal and background count 
distributions for all the dSph galaxies that were analyzed. 
The background counts have been scaled taking into 
account the solid angle ratio between the signal and 
background regions, according to Eq. [5] (Eq. [TT] in the 
case of the stacking analysis). In all cases no evidence 
is observed of a net signal excess with respect to the 
background in any energy bin. 

Fig. [4] shows the upper limits at 95% CL on the flux 
(top panel) and on & PP (E) (bottom panel) as function 
of the energy, for each of the dSph galaxies considered 
in this analysis. The more constraining limits are those 
obtained from the dSph galaxies with the highest J- 
factors. 

Once the individual sources were analyzed, the 
stacking analysis and the composite analysis were 
implemented following the procedures described in ijIIIBI 
and in mil CI As mentioned above, for the stacking 
analysis the counts from the signal and background 
regions corresponding to each source were added, and the 
upper limits on the signal were evaluated following the 
same procedure as for individual sources. This approach 
is equivalent to stacking the events from all the dSph 
galaxies and then analyzing the single image obtained 
from the superposition of all the individual images (see 
the last plots in Figs. [2] and [3]). 

It is worth noting that in the high-energy bands, where 
the counts in the signal and background regions are both 
null (i.e. n = m = 0), as shown in Fig. [3l the upper 
limit evaluated on the signal counts is always constant 
(i.e., at 95% CL the upper limits on s corresponds to 
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FIG. 2: Photon count maps in the observed energy range from 562 MeV to 562 GeV for the dSph galaxies considered in this 
analysis. The black circles indicate the cones of angular radii of 0.5°, 5° and 6°, representing the boundaries of the signal and 
background regions. The sources in the 2FGL Catalog are indicated with crosses. Each map is centered on the position of the 
corresponding source. The map in the bottom right panel was obtained by stacking the data from all the dSph galaxies. 
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FIG. 3: Count distributions in the observed energy range from 562 MeV to 562 GeV for the dSph galaxies considered in 
this analysis. The black points represent the counts in the signal region; the grey areas represent the equivalent number of 
background counts (i.e., the counts in the background region scaled by the coefficient c) with their errors. The bottom right 
panel shows the distribution obtained by stacking all the dSph galaxies. 
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about 3 counts). Hence, when performing the stacking 
analysis, the upper limit on the flux (see top panel of 
Fig. H|) will decrease linearly with the number of stacked 
sources. In other words, in the stacking analysis the 
observations of different sources are added and the result 
is expected to be equivalent to a single observation of an 
individual source with a total live time corresponding to 
the sum of the live times of each observation. On the 
other hand, in the low-energy band, where the counts 
in the signal region are roughly equal to those in the 
background region (i.e. n ~ c to), the upper limit on the 
signal counts s is roughly proportional to the square root 
of the observed events [15|]. In this case, since the total 
live time will be roughly proportional to the number of 
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FIG. 4: Top panel: Upper limits at 95% CL on the gamma-ray 
flux as function of energy. Bottom panel: Upper limits at 95% 
CL on & PP (E) as function of energy. The colored symbols 
correspond to the results obtained from the individual dSph 
galaxies. The open black circles indicate the results obtained 
from the composite analysis, while the filled black circles 
indicate the results obtained from the stacking analysis of 
all the dSph galaxies. 



stacked sources, the upper limit on the stacked flux will 
improve with the square root of the number of stacked 
sources. 

In the bottom panel of Fig. [4] the upper limits on 
$pp(E) are shown for all the candidate sources as well 
as for the stacking and composite analyses. As shown 
in Tab. U the J-factors of the 10 dSph galaxies studied 
in the present analysis are distributed in an interval 
that spans two orders of magnitude. As a consequence, 
since the upper limits on the photon fluxes are roughly 
similar, the upper limits on $> PP (E) will span two 
orders of magnitude. The stacking analysis improves 
the upper limits on & PP (E) by a factor that ranges 
from a few to about 10 with respect to those obtained 
from the analysis of individual dSph galaxies. When 
considering the quantity § PP (E), since the J- factor used 
in the stacking analysis is evaluated as the average of 
the J-factors of individual sources weighted with their 
exposures, the result is an improvement of a factor of a 
few with respect to the upper limits obtained from the 
analysis of the source with the highest J-factor. 

The results from the composite analysis are in general 
more constraining than those from the stacking analysis. 
This could be due to the fact that the random variable 
used to evaluate the upper limits in the composite 
analysis is oc \\ i Si/Ji 1 while the random variable used 
in the stacking analysis is oc Q^ Si)/ (J). This means 
that the "effective J-factor" for the composite analysis 
could in principle be different from that for the stacking 
analysis. 

The measured upper limits on $> pp (E) were converted 
into upper limits on (av) following the procedure 
described in §111 Al Fig. [5] shows an example of this 
calculation in the case of Segue I for the annihilation 
channels ^ + /i~, t + t~, bb and W + W~ . The upper 
limits were also evaluated taking into account separately 
the uncertainties on the effective area and on the J- 
factors. To describe the systematic uncertainties on the 
effective area we assumed a uniform PDF centered on the 
average value A(E) in each energy bin with fluctuations 
of ±10%. These uncertainties have a negligible effect on 
the upper limits. As discussed above, the effects on the 
upper limits due to the systematic uncertainties on the 
J-factor were evaluated assuming for J a uniform PDF 
in a range corresponding to the 68% area of the actual 
J-factor distribution. As shown in Fig. [SJ the effects of 
the uncertainties on the J-factor are not negligible and, 
depending on the source under investigation, the upper 
limits on § PP (E) can increase by up to a factor of a few. 

Fig. |6] shows the upper limits at 95% CL on (av) 
obtained from the analysis of individual dSph galaxies, 
from the stacking analysis and from the composite 
analysis, as a function of the WIMP mass for the 
annihilation channels n + n~ , r + r _ , bb and W + W~ . The 
upper limits obtained by taking into account the effects 
of the uncertainties on the J-factors are also shown. 
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FIG. 5: Evaluation of the upper limits on (av) as a function of the true energy for several WIMP mass values in the case 
of Segue I. The blacks lines with full circles correspond to the upper limits at 95% CL on & PP (E). The colored lines, each 
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upper limits. The four panels refer to WIMP annihilation into fi + /j,~ , t + t~ , bb and W + W~ (as labeled). The dashed lines with 
filled squares and the dotted lines with filled triangles indicate the upper limits evaluated taking into account the systematic 
uncertainties on the effective area and on the J-factors, respectively. The effects of the systematic uncertainties on the effective 
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V. ANALYSIS OF THE MILKY WAY HALO 

The study of the Milky Way halo is quite complex 
because its gamma-ray emission has to be disentangled 
from that of known gamma-ray sources. However, a 
possible approach to the study of the Milky Way halo 
is that of selecting a set of sky directions that are well- 
separated from known gamma-ray sources. 

For this analysis a set of 1000 random directions 
was generated in the sky, each direction located at an 
angular distance of at least 3° from all the 1873 point 
sources and at least 3° plus twice the size of the semi- 
major axis from all the 11 extended sources in the 2FGL 
Catalog [la ]. The random positions are illustrated in 
Fig. [7] in Galactic coordinates. Since many gamma-ray 
sources are concentrated in the region of the Galactic 
plane, we also decided to perform a separate analysis 
selecting only random directions at an angular distance 



larger than 10° from the Galactic plane (i.e. all directions 
with Galactic latitude \f3\ > 10°). A subset of the initial 
sample, consisting of 866 random directions, was used for 
this analysis. 

The analysis of the Milky Way Halo was performed 
selecting P7CLEAN_V6 class events in order to guarantee 
optimal rejection of the charged particle background. 
As in the case of the dSph galaxies, discussed in mVl 
the data analysis was performed selecting gamma rays 
with energies from 562 MeV to 562 GeV, with the energy 
interval being divided into 12 bins, equally spaced on a 
logarithmic scale. 

In the case of the Galactic halo analysis, we 
hypothesized an extreme scenario in which all the 
detected photons originate from DM annihilation. In this 
analysis the upper limits on $> PP (E), and consequently 
those on (av), were evaluated assuming an absence of 
background events. For each random direction the signal 
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FIG. 6: Upper limits at 95% CL on {av) as function of the WIMP mass for the annihilation channels fi + \i~ , t + t~ , bb and 
W + W~ . The plots show the upper limits obtained from the analysis of individual dSph galaxies and from the stacking and 
composite analyses. The continuous lines indicate the upper limits obtained neglecting the systematic uncertainties, while the 
dotted lines indicate the upper limits obtained including the uncertainties on the J-factors. The long dashed line corresponds 
to the canonical value of the annihilation cross section of 3 x 10 -26 cm 3 s _1 in the thermal relic WIMP scenario. 



region was defined as a cone of angular radius A8 = 
1° centered on it. The analysis was then performed 
stacking the data from all the random directions, without 
background subtraction (i.e. c = and m = 0). This 
strategy relies on the assumption that a possible DM- 
induced gamma-ray flux must be lower than the total 
observed flux. The upper limits obtained from this 
analysis therefore will be conservative. 

For any given random direction (hereafter we will 
refer to random directions as sources) the J-factor was 



evaluated using Eq. [2] Since the signal region is a narrow 
cone of 1° angular radius, Eq. [5] reduces to: 



J k AVI / p 2 (l{4>))dl. 



(34) 



where Ail w 9.6 • 10~ 4 sr is the solid angle corresponding 
to the signal region. In the previous equation we 
explicitly wrote the dependence of the DM density on 
the angle if), which represents the angular separation of 
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FIG. 7: Count map used for the analysis of the Milky Way 
Halo. The map was built in the Galactic reference frame, 
using the HEALPix pixelization scheme with N s id e — 128 
(196608 pixels, each covering a solid angle of 6.4 ■ 10~ 5 sr), 
and is displayed in the Aitoff projection. The Galactic Center 
is in the middle of the map. The 1000 random directions 
are indicated with the red markers; the sources in the 2FGL 
Catalog are indicated with the blue markers. 
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For the DM density we assumed a Navarro-Frenk- 
White (NFW) profile [22J: 



p(r) 



Pa 



(r/r s )(l + r/r s y 



(36) 



where po = 0.3GeV/cm 3 and r s — 20kpc. The 
coordinate r in Eq. (35] represents the distance from the 
Galactic Center, and is given by: 



l 2 + 



R 2 



2lR cost(j 



(37) 



where Rq — 8.5 kpc is the distance of the Galactic Center 
from the Sun and / is the distance of the line element dl 
from the Sun. 

Fig. [5] shows the distribution of the J- factors evaluated 
for the sky directions shown in Fig. [7] As expected, the 
directions with the highest J-factors are near the Galactic 
Center; on the other hand, the directions with the lowest 
J-factors are near the Galactic Anti-center. 

The stacking analysis was performed as described 
in §111 Bl In our model we assumed that all photons 
originate from DM, and the upper limits on the signal 
counts were evaluated using the PDF in Eq. [5J We also 
performed the composite analysis of the random sources 
following the procedure described in £|IIIC[ and in this 
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FIG. 8: Distribution of the J-factors evaluated for the random 
directions in the Galactic halo. The dashed filled region 
shows the J-factors of all the 1000 directions; the grey filled 
region shows the J-factors for the 866 sky directions that are 
separated more than 10° from the Galactic plane. 



the source from the Galactic Center. 

In performing the calculations we assumed that the 
detector is located at the Sun's position and we used 
the Galactic reference frame. The angle ip can then be 
calculated from the Galactic longitude and latitude (A, 
0) using the following relation: 



cos?/' = cos A cos p. 



(35) 
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FIG. 9: Comparison of the 95% CL upper limits on $> PP (E) 
as function of the energy for all of the 1000 random directions 
and for the 866 directions outside the Galactic plane (|/3| > 
10°). The upper limits obtained from the directions 
corresponding to the highest and lowest J-factors are also 
shown. The lowest J- factor among the 1000 sky directions 
is associated with the direction (in Galactic coordinates) 
(A,/3) = (183.616°, -0.189°), while the highest J-factor is 
associated with the direction (356.796°, 9.23°). The lowest 
J-factor among the 866 sky directions outside the Galactic 
plane is associated with the direction (179.968°, —11.5651°), 
while the highest J-factor is associated with the direction 
(0.621474°, 10.0793°). 
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FIG. 10: Evaluation ol the upper limits on (crv) as a function of the true energy for several WIMP mass values from the 
Galactic halo analysis with the stacking of all the 1000 random directions. The black lines correspond to the upper limits at 
95% CL on <3> (E). The colored lines, each corresponding to a given value of the WIMP mass, indicate the maximum allowed 
values of & PP (E) that do not exceed the measured upper limits. The four panels refer to the WIMP annihilations into fi + fj,~ , 
t + t~ , bb and W + W~ , respectively, as labeled. 
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FIG. 11: Upper limits at 95% CL on (av) as a function of the WIMP mass for the annihilation channels fi + fi~ , t + t~, bb and 
W + W~ . The plot shows the results obtained from the analyses of the Milky Way halo with all the 1000 directions (left panel) 
and with the 866 directions with \/3\ > 10° (right panel). The dashed line is the annihilation cross section of 3 x 10~ 
in the canonical thermal relic WIMP scenario. 
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case the results were equivalent to those obtained from 
the stacking analysis. 

Fig. H] shows the upper limits at 95% CL on $ pp (E) 
for the Milky Way halo evaluated using all the 1000 
random directions and only the 866 directions with \j3\ > 
10°. The upper limits obtained from the directions with 
the lowest and the highest J-factors are also shown. 
The directions with the highest J-factor yield the more 
constraining upper limits on $ PP (_E). This result is 
not completely obvious, since the directions with the 
higher J-factors are in the region of the Galactic Center, 
where a high number of photons is expected, while the 
directions with the lower J-factors are in the region of the 
Galactic Anti-center, where a lower number of photons 
is expected. The upper limits obtained from the analysis 
of the direction with the highest J-factor are also more 
constraining than the ones from the stacking analysis, 
with the exception of the high energy regime, where the 
constraints from the stacking analysis are tighter. This is 
due to the fact that in the high-energy regime the number 
of events in the signal regions is low, and consequently 
the upper limits on the fluxes decrease with increasing 
live time. 

Fig. [10] shows the procedure used to convert the 
measured upper limits on $ PP (_E) into upper limits on 
(cry) in the case of the stacking analysis of all the 1000 
random sources. As in the case of the dSph analysis, we 
imposed the requirement that the flux values predicted 
from the DM annihilation scenarios must not exceed the 
measured upper limits in any energy bin. 

In Fig. |TT| the upper limits at 95% CL on (av) as 
a function of the WIMP mass are shown for WIMP 
annihilations into /x + /i~, t + t~, bb and W + W~ for the 
Milky Way halo. These results have been obtained from 
the analysis of all the 1000 random directions (left plot) 
and from the analysis of the 866 directions with \j3\ > 10° 
(right plot) (see Fig. |H]). If only the direction with 
the highest J-factor were analyzed, the upper limits on 
(av) in the low WIMP-mass regime would be about a 
factor 10 more constraining than those obtained from 
the combined analyses. 



VI. DISCUSSION 

In this work we used the gamma-ray data collected 
by the Fermi LAT during its first 3 years of operation 
to set constraints on the parameter (av) assuming DM 
annihilation into various channels. 

We studied a set of 10 Milky Way dSph satellite 
galaxies and the Milky Way halo. The dSph 

galaxies were analyzed both individually and collectively, 
implementing dedicated stacking and composite analysis 
procedures. The Milky Way halo was studied by 
randomly sampling a set of 1000 sky directions well- 
separated from all the gamma-ray sources of the 2FGL 
Catalog and performing a stacking analysis. The 
data analysis was performed using a model-independent 



method that allows upper limits to be set on the gamma- 
ray fluxes starting from the observed events using a 
Bayesian approach. The constraints on (av) were derived 
requiring that the predicted fluxes from the models must 
not exceed the measured ones. 

The analysis of the dSph galaxies yields upper limits 
on (av) that are lower with respect to the predictions 
from a canonical thermal WIMP scenario for the t + t~ 
and bb final states up to masses of few tens of GeV. 
This is found in the stacking and in the composite 
analysis results, but also in the results of the individual 
analyses of the dSph galaxies with the highest J-factors. 
The uncertainties in the J-factor calculation and on 
the effective area of the LAT were also included in the 
present analysis, and do not affect significantly the upper 
limits. Our results are consistent with recent analyses 
yl Q performed using different approaches. However, 
we emphasize that the upper limits on the parameter 
(av) depend strongly on the values of the J-factor. In 
particular, since no evidence of a gamma-ray flux is 
observed from any dSph galaxy, the upper limits on 
<& PP (E) and consequently those on (av) will scale with 
the J-factor (see Fig. [4]). 

For comparison, the analysis of the dSph galaxies was 
also performed using the P6_V3_DIFFUSE IRFs, as were 
used in Ref. j5j, and the results were found to be in 
agreement with the ones already presented here, which 
were obtained using the P7SOURCE_V6 IRFs. 

The analysis of the Milky Way halo, performed looking 
at a set of 1000 clean sky directions, yields upper limits on 
(av) that range from 10~ 25 to 10~ 23 cm 3 jr 1 for WIMP 
masses below lOTeV or more for the bb and W + W~ 
annihilation channels. More constraining limits on (av) 
can be obtained in the low-energy region if the analysis is 
limited to the direction with the highest J-factor. These 
limits were evaluated assuming a NFW profile with a 
DM density at the solar circle of 0.3 GeVcm -3 , and their 
values depend on the DM mass density profile. A recent 
analysis suggested a revised value of the DM density 
at the solar circle of 0.43GeVcm -3 [23J]; under this 
assumption the upper limits on (av) would improve by a 
factor of 2. Nevertheless we note that the sky directions 
used for the present analysis are sufficiently far away from 
the Galactic Center that the J-factors evaluated with 
different DM density profiles would likely yield similar 
results (see Ref. [24j ). The current results are consistent 
with the results obtained by Ref. [25|, where an analysis 
of the all-sky Fermi LAT data was performed with a 
different model-independent technique. The Milky Way 
halo can also be studied following a different approach, in 
which a model is assumed for the Galactic diffuse gamma- 
ray emission and a fit of the DM signal together with 
the diffuse component is performed (e.g., see 26]). Our 
results are also consistent with those obtained from that 
analysis. 

Both from the limits on the dSph galaxies (Fig. [S]) and 
those on the Milky Way halo (Fig. (TTJ, it is possible 
to restrict the range of allowed WIMP masses assuming 
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the standard thermal relic scenario. We also note that 
the DMFIT package may underestimate the gamma- 
ray fluxes for WIMP masses above 1 TeV, since it does 
not include radiative electroweak corrections; hence, the 
limits on (av) for large masses can be viewed as more 
conservative than those in the low mass region. 



VII. CONCLUSIONS 

We developed a model-independent approach to set 
upper limits on the energy spectra of both individual 
and multiple gamma-ray sources using the data collected 
by the Fermi LAT. In this paper we presented the 
results obtained from the application of this technique 
to the study of a set of dSph galaxies and to the study 
of the Milky Way halo. These results were used to 
derive constraints on DM annihilation cross sections 
into different channels. We emphasize that the analysis 
techniques presented in this paper are general, and are 
suitable for applications where the study of a class of 
sources, even faint sources, with common features have 
to be studied. 

The data analysis technique illustrated in the present 
paper allows us to set robust upper limits on the energy 
spectra of candidate gamma-ray sources. The upper 
limits on the gamma-ray fluxes are in fact derived 
starting from the data, without assuming any model 
for the background and for the source spectral shapes. 
The signal is evaluated selecting events from a cone 
centered on the source position, while the background 
is evaluated selecting events in a region close to the 
source under investigation. Both signal and background 
fluctuations are described in the framework of Poisson 
statistics, and the upper limits on the signal counts, and 
consequently on the flux, are computed following the 
Bayesian approach. A stacking analysis and a composite 
analysis procedure have also been developed to perform 
the collective study of multiple candidate sources with 
common features. 

The analysis methods presented in this paper can 
also be applied when several measurements of a given 
physical quantity, each one resulting into a confidence 
interval, have to be combined into a unique confidence 
interval taking all the results into account. The 
systematic uncertainties can also be incorporated in the 
analysis by introducing proper nuisance parameters in 
the probability distribution functions. 
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Appendix A: Derivation of Equation 1121 

When performing the stacking analysis of a set of 
multiple sources, the counts from the individual signal 
and background regions are added. The total counts in 
the signal and background region are then given by: 



n = > n 



E 

i 

E 



(Al) 
(A2) 



where rii and rrii are respectively the counts in the signal 
and in the background region of the «-th source. 

According to the assumptions in ^III At n, and m, 
are both Poissonian with expectation values Si + Cibi 
and bi (in the following we shall use the notation n% ~ 
V(si + Cibi) and rrii ~ V(bi)). Hence, from the definitions 
of n and m it follows that: 



(A3) 

(A4) 



In the stacking analysis the true values of the signal 
and of the background counts are defined as: 




s = E s > 

i 

& = $> 



(A5) 
(A6) 
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This definition automatically implies that m ~ V(b). 
However, to ensure that n ~ V(s + cb), the coefficient c 
must be defined as: 



The values of bi are not known, but they can be replaced 
with their best estimators b* = rrii + 1 [23]. In this way, 
Eg. 1X71 reduces to Eg. [HI 






(A7) 
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